/*Program for estimating betak and betay using LP variables (capital, materials 
	and industry output)*/
	
/*First stage*/
*Create polynomial of degree `pmax' in k, y, and m
local i = 1
	foreach x in k m q {
		quietly gen double var_`i' = `x'
		local i = `i'+1
		}

	forv i = 1/3 {
		forv j = `i'/3 {
			quietly gen double var_`i'`j' = var_`i'*var_`j'
			}
		}
	
*Run first-stage regression
qui reg y_rev_net var*

*Predict phihat
qui predict phihat
drop var*

/*Second stage*/
*GMM program
capture program drop gmm_ky_LP_con
	program define gmm_ky_LP_con
	tempvar e omega omlag omlag2 omlag3 epsilon mom1 mom2 xx
	matrix score `e'=`1'
	qui gen double `omega' = phihat-`e'
	sort lbd year
	qui gen double `omlag'=L.`omega'
	qui gen double `omlag2'=`omlag'^2
	qui gen double `omlag3'=`omlag'^3

	qui reg `omega' `omlag' `omlag2' `omlag3'
	qui predict `xx'
	qui gen `epsilon'=`omega'-`xx'

		qui gen double `mom1' = (`epsilon'*k)
		qui su `mom1'
		scalar sumom1 = (r(sum))^2
		
		qui gen double `mom2' = (`epsilon'*q)
		qui su `mom2'
		scalar sumom2 = (r(sum))^2
		
		scalar c1 = `1'[1,1]+$aLprog
		scalar c2 = `1'[1,1]
		scalar c3 = `1'[1,2]
		
		if c1 > .99 {
			scalar obj = -999999999999
			}
		else if c2 < 0.01 {
			scalar obj = -999999999999
			}
		else if c3 < 0.01 {
			scalar obj = -999999999999
			}
		else if c3 > 0.5 {
			scalar obj = -999999999999
			}
		else {
			scalar obj = -(sumom1 + sumom2)
			}
		scalar `2'=obj
		
	end
	
*Initial guess
reghdfe y_rev_net k q, ab(lbd)
mat initols = e(b)
	/*Enforcing that initial conditions meet the constraint criteria*/
	local akmat = initols[1,1]
	if `akmat' < 0.01 {
		mat initols[1,1] = 0.5
		}
				
	local summat = $aLprog + initols[1,1]
	if `summat' > .99 {
		mat initols[1,1] = 0.98-$aLprog
		}
			
	local qmat = initols[1,2]
	if `qmat' < 0.01 {
		mat initols[1,2] = 0.25
		}
	if `qmat' > 0.5 {
		mat initols[1,2] = 0.25
		}
		
*Optimization
amoeba gmm_ky_LP_con initols obj vcoef . 200 0.0000001
